{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parks–McClellan filter design algorithm Book:\n",
    "\n",
    "  https://www.detailedpedia.com/wiki-Remez_exchange_algorithm\n",
    "  \n",
    "  https://www.detailedpedia.com/wiki-Parks%E2%80%93McClellan_filter_design_algorithm\n",
    "\n",
    "Digital filter:\n",
    "\n",
    "  https://www.detailedpedia.com/wiki-Digital_filter\n",
    "\n",
    "Butterworth filter (detailedpedia.com):\n",
    "\n",
    "  https://www.detailedpedia.com/wiki-Butterworth_filter\n",
    "\n",
    "Remez from scipy:\n",
    "    https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.remez.html\n",
    "\n",
    "斯坦福 Filter 教程 by Jos:\n",
    "https://ccrma.stanford.edu/~jos/filters/\n",
    "\n",
    "DADiSP Online Help (SOS2TF):\n",
    "https://www.dadisp.com/webhelp/mergedProjects/refman2/FncrefSZ/SOS2TF.htm\n",
    "\n",
    "THE FOUR DIRECT FORMS: https://www.dsprelated.com/freebooks/filters/Four_Direct_Forms.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 22050   # Sample rate, Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_response(w, h, title):\n",
    "    \"Utility function to plot response functions\"\n",
    "    fig = plt.figure(dpi=800)\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.plot(w, 20*np.log10(np.abs(h)))\n",
    "    ax.set_ylim(-40, 5)\n",
    "    ax.grid(True)\n",
    "    ax.set_xlabel('Frequency (Hz)')\n",
    "    ax.set_ylabel('Gain (dB)')\n",
    "    ax.set_title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first example is a low-pass filter, with cutoff frequency 8 kHz. The filter length is 325, and the transition width from pass to stop is 100 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff = 8000.0    # Desired cutoff frequency, Hz\n",
    "trans_width = 100  # Width of transition from pass to stop, Hz\n",
    "numtaps = 325      # Size of the FIR filter.\n",
    "taps = signal.remez(numtaps, [0, cutoff, cutoff + trans_width, 0.5*fs],\n",
    "                    [1, 0], fs=fs)\n",
    "w, h = signal.freqz(taps, [1], worN=2000, fs=fs)\n",
    "#plt.figure(1, dpi=800)\n",
    "plot_response(w, h, \"Low-pass Filter\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows a high-pass filter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff = 2000.0    # Desired cutoff frequency, Hz\n",
    "trans_width = 250  # Width of transition from pass to stop, Hz\n",
    "numtaps = 125      # Size of the FIR filter.\n",
    "taps = signal.remez(numtaps, [0, cutoff - trans_width, cutoff, 0.5*fs],\n",
    "                    [0, 1], fs=fs)\n",
    "w, h = signal.freqz(taps, [1], worN=2000, fs=fs)\n",
    "plot_response(w, h, \"High-pass Filter\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows a band-pass filter with a pass-band from 2 kHz to 5 kHz. The transition width is 260 Hz and the length of the filter is 63, which is smaller than in the other examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band = [2000, 5000]  # Desired pass band, Hz\n",
    "trans_width = 260    # Width of transition from pass to stop, Hz\n",
    "numtaps = 63         # Size of the FIR filter.\n",
    "edges = [0, band[0] - trans_width, band[0], band[1],\n",
    "         band[1] + trans_width, 0.5*fs]\n",
    "taps = signal.remez(numtaps, edges, [0, 1, 0], fs=fs)\n",
    "w, h = signal.freqz(taps, [1], worN=2000, fs=fs)\n",
    "plot_response(w, h, \"Band-pass Filter\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The low order leads to higher ripple and less steep transitions.\n",
    "\n",
    "The next example shows a band-stop filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band = [6000, 8000]  # Desired stop band, Hz\n",
    "trans_width = 200    # Width of transition from pass to stop, Hz\n",
    "numtaps = 175        # Size of the FIR filter.\n",
    "edges = [0, band[0] - trans_width, band[0], band[1],\n",
    "         band[1] + trans_width, 0.5*fs]\n",
    "taps = signal.remez(numtaps, edges, [1, 0, 1], fs=fs)\n",
    "w, h = signal.freqz(taps, [1], worN=2000, fs=fs)\n",
    "plot_response(w, h, \"Band-stop Filter\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import sos2tf\n",
    "from scipy import signal\n",
    "\n",
    "sos1 = signal.butter(2, 0.1, btype='low', output='sos')\n",
    "print('sos1 = ', sos1)\n",
    "b1, a1 = sos2tf(sos1)\n",
    "print(\"b1, a1 =\", (b1, a1))\n",
    "\n",
    "sos2 = signal.butter(2, [0.04,0.1], btype='band', output='sos')\n",
    "print('sos2 = ', sos2)\n",
    "b2, a2 = sos2tf(sos2)\n",
    "print(\"b2, a2 =\", (b2, a2))\n",
    "\n",
    "# sos4 = signal.butter(4, [0.04,0.1], btype='band', output='sos')\n",
    "# print('sos4 = ', sos4)\n",
    "# b4, a4 = sos2tf(sos4)\n",
    "# print(\"b2, a2 =\", (b2, a2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.polynomial import Polynomial\n",
    "np.polynomial.set_default_printstyle('unicode')\n",
    "\n",
    "# B1 = Polynomial([ 0.00782021,  0.01564042,  0.00782021])\n",
    "# B2 = Polynomial([ 1.,         -2.,          1])\n",
    "# A1 = Polynomial([1.,         -1.76530323,  0.83753376])\n",
    "# A2 = Polynomial([1.,         -1.89572603,  0.91459788])\n",
    "B1 = Polynomial(sos2[0][:3])\n",
    "B2 = Polynomial(sos2[1][:3])\n",
    "A1 = Polynomial(sos2[0][-3:])\n",
    "A2 = Polynomial(sos2[1][-3:])\n",
    "print(\"(B1, A1) = \", (B1, A1));\n",
    "print(\"(B2, A2) = \", (B2, A2));\n",
    "B = (B1*B2).coef\n",
    "A = (A1*A2).coef\n",
    "print('B, A = ', (B, A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal\n",
    "#from scipy.signal import butter\n",
    "order=2\n",
    "sos1 = signal.butter(order, 0.1, btype='low', output='sos');\n",
    "sos2 = signal.butter(order, [0.04,0.1], btype='band', output='sos');\n",
    "b, a = signal.sos2tf(sos2)\n",
    "print('b, a', (b, a))\n",
    "z, p, k = signal.tf2zpk(b, a)\n",
    "print('z, p, k = ', (z, p, k))\n",
    "sos3 = signal.butter(order, 0.04, btype='high', output='sos')\n",
    "print('sos1: ', sos1)\n",
    "print('sos2: ', sos2)\n",
    "print('sos3: ', sos3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('len(sos2)', len(sos2))\n",
    "print('sos2.shape', sos2.shape)\n",
    "print('sos2[:, 2]', sos2[:, 2])\n",
    "print('sos2[:, 5]', sos2[:, 5])\n",
    "print('(sos2[:, 2]==0)',(sos2[:, 2]==0).sum())\n",
    "print('(sos2[:, 5]==0)',(sos2[:, 5]==0).sum())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
